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I.  Introduction 


The  purpose  of  this  project  was  to  determine  false  alrarm  rates  and  pro¬ 
bability  of  detection  for  the  use  of  Fourier  descriptors  in  recognition  of  tar 
gets  in  aerial  imagery.  Another  goal  was  to  develop  the  software  package 
for  Fourier  descriptor  analysis  so  that  it  could  be  installed  on  a  VAXL1/780 
computer  ruining  the  UNIX  operating  system  at  Feme  Air  Development  Center. 


II.  Fourier  Descriptor  Theory 

The  Fourier  descriptor  is  a  method  using  basis  functions  to  decompose 
or  characterize  the  shape  of  an  object.  The  basis  functions  of  the  Fourier 
series  have  nice  properties .  All  the  rigid  motions  in  the  plane  and  dilations 
can  be  performed  by  simple  operations  in  the  representation  space.  By 
tracing  the  boundary  of  the  silhouette,  a  complex  function  of  time  is 
obtained.  This  function  is  periodic  and  can  be  expanded  in  a  Fourier  series. 
Several  authors  have  extended  the  list  of  known  properties  and  methods  of 
normalization  and  recognition  [GRAN72],  [RICH74],  [PERS77],  [WALL80]. 

Given  the  silhouette  of  the  object,  trace  the  boundary  counter¬ 
clockwise  at  a  uniform  velocity.  After  one  complete  circuit  around  the 
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boundary  the  function  repeats.  Let, 

7(0  &  (*<0  .  v(0)  =  *(0  +  <  v(Oc«xR  . 

v(<)  i  velocity  * 

ot 

So, 

l„(OI**  MIL ?.7*(0 .. 

1  '  Ot  ot 


A  uniform  tracing  velocity  implies  that 

|v(t)|  ®  V*(<)z  +  fjr(f)*  =  Constant. 

Since  y(t)  is  periodic 

7  <t)=7(t  +  n. 


where  T  =  period  and  vT  -  L  -  arc-length.  For  simplicity  let  v  =  1  space 
unit  /  time  unit.  Since  y(t)  is  periodic,  we  can  expand  in  its  Fourier  series 


2wn , 


r(0  =  E  c»  « 

»■  — «■ 


where 


.SlHlj 


iq+r 

C ,  =  V/l’t')*"  r  «“• 


All  the  well  known  properties  of  Fourier  series  analysis  are  now  at  our 
disposal  in  investigating  the  properties  of  this  representation. 

^  2*71 1 

Convergence:  First,  the  basis  functions  e  7  forms  an  orthonormal 
basis  complete  in  Lz.  So,  for  any  curve  that  is  connected  and 
rectifiable',  the  series  converges. 

Continuity:  y(t)  is  everywhere  continuous.  So,  cn  ~ 


D 


2) 

3) 


n‘ 


Translation: 


y(0  =  7<0+Z(r*  c'0  =  c 0+Z0 


4)  Rotation :  Rotation  of  the  shape  through  an  angle  a 


Y(t)  =  7 (0«“  -  c'»  =  e»  •**. 


5)  Dilation  (Seals): 

Y(t)  =  *7(f)  -*  c'n  =Ac„. 


6)  Starting  point:  For  a  shift  in  starting  point, 
y(0  =  jit ~ *o)  -*  c'n  =  cn  e  T  ®. 


7)  N-fold  symmetry:  For  a  shape  having  N-fold  rotational  symmetry, 

<2SL  _  iSZL  .awn  T 

y(t)s  N  —  y(t  +  jj)  -  cn«  N  =  cne  T  N  . 

We  have 

Cn(l“«  "  )  =0. 


So,  for  each  n  either  cn  =  0  or 
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This  means  that  (n-l)  must  be  a  multiple  of  N. 

Hence,  for  an  object  to  have  N-fold  rotational  symmetry,  the  only  non¬ 
zero  coefficients,  cm.  are  those  such  that  n  =  pN+ 1,  for  some  integer  p . 

8)  Axial  symmetry:  For  a  shape  which  is  symmetric  through  its  centroid, 
y(tt)  =  -7«+f*  • 


So, 


Now, 

(1  -  e"*1*)  =  2,  n  odd. 


=  0,  n  even. 


Hence, 


*0'  2  2m 

cn-yf  7(0«”*  T  *dt  ,n  odd. 

1  <o 

=  0  ,  n  even. 

Thus,  a  shape  exhibiting  axial  symmetry  has  only  odd  non-zero 
coefficients. 


9)  Reflections:  Let  y'  be  the  reflection  of  y  about  the  vertical  axis  through 
the  centroid.  Hence,  we  have 

y(t)  =  -*<-<  +  +  iy(-t  + 

Writing 


iSSX-t 

*(0  =  2«n«  T 


iSELi 

y(t)  =  2X*  T  • 


From  our  previous  definition  of  cn, 

Cn  =  <*n  +  *•  =  One4(n+1>w  +  ^*0**”. 


This  gives  us  the  following  relationship  between  y  and  y': 
k'_*  |  =  |cn|  and£c  =  (n  +  l)rr  -  Lcn  . 

11)  i4reo.-  For  a  simple  closed  curve  (not  necessarily  traced  at  an  uniform 
velocity), 

Area  ~  xdy  -  )-f  ydx  . 

C  7  C  7 

<0+rf  1 

Area  =  4“/  “  *y|  dt 

c  <0 

Area  =  -r  E  ™  |cn  |2 


When  ths  curve  is  traced  at  a  constant  speed. 


Ana  s  jm  |cn|*- 

12)  Line  shop*:  For  a  line  shape,  i.e.  one  with  no  interior  7(0  =  y(T-t). 
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Let  X=7W  in  second  integral  and  substitute  7(X)  for  7(f -X)  to  obtain 
T_ 

cn  =  j-fy(t)cos(^p-t)dt  ->  cn  =  -  Area  *  0. 

-  An  equivalent  Fourier  series  representation  of  7(0  is  the  following 
[KUHL] : 

t  1 

*(0  =  2  Kco*(2ES.O 

n»0l  1  1 

V (0  =  2  Z^rtn(2ES-0 

»»o[  *  1 

Now,  y{t)  =  x(0  +  iy(t).  So,  substituting  the  above  for  *  and  y 

r(0  =  )J)jke“<SrL*) +  Jt-M4?*-*) 

+  <jc;co*(SEt.o  +  n,«*»(*5-jr  ■ 


c 


=  j{(4.  -  A.)  +  i(B»  +  CW)]  ,  n=  1.2,  •  •  ■ 


The  Normalized  Fourier  Descriptor 

In  this  section  we  will  discuss  the  procedure  used  to  normalize  the 
characterization  of  the  shape  to  a  standard  orientation,  size,  and  starting 
point.  To  compare  two  Fourier  descriptors  (FD’s)  for  shape  recognition  we 
need  to  scale,  rotate,  and  shift  the  starting  point  of  the  unknown  object 
boundary  curve  to  that  of  the  template  curve. 

Persooi?  &  Fu  [PERS77]  describe  a  method  that  is  "optimum"  in  that  it 
minimizes  the  mean-squared  error  between  the  unknown  and  the  template 
shape.  Given  two  shapes  y  and  y'  having  coefficients  cn  and  c'n,  respec¬ 
tively,  minimize  the  distance 


d(y.Y)  = 


1 

2 


In  practice,  only  N  coefficients  of  the  FD  for  the  shape  are  retained.  We 
want  to  scale,  rotate,  and  shift  the  starting  pint  by  introducing  the  parame¬ 
ters  s ,  a,  tp  ,  respectively.  Hence,  we  want  to  minimize  the  quantity 

n*—N/2,  n#0L  1 


Letting 

cnc  '»  =Pnelf" 

we  have 

£c»cn*  +  s2£c'nc',*  -  2sZpncos(1rn+n<p+<x)  . 

To  minimize  this  expression  with  respect  to  s,  a,  ^  take  partial  deriva¬ 
tives  and  set  them  equal  to  zero. 

-£-=  2s2c'nc',f  -  2 Zpncos(itn+n<p+a) 


-9- 


^-=  2s^pnsin(i/n+n(p+a) 


“-=  2 sSpnnsin^+np+a)  . 


By  setting  the  partials  equal  to  zero  the  following  equations  are  obtained 


%pncos  ( j/n+ntp+<x ) 
£c'«c'n* 


(1) 


tan(<p)  =  — 


£pnsm(y>n 

SpnCOsCV'n+ny) 


(2) 


ian(^)  =  - 


%pwn  sm(yn->-ny) 
pnn  cos(^n+^V) 


(3) 


Combining  (1)  and  (2)  into  one  equation  in  <p 

f(<p)  =  £pnsin(tn+7l^)S7lPncos(V,n+n«») 


-  Zpncos(yn+n<p)npnsintyn+n(p)  . 


The  optimum  value  for  <p  is  obtained  when  /  (<p)  is  zero.  The  corresponding 
values  of  s  and  a  can  be  obtained  from  equations  (1)  and  (2). 

Because  only  N  coefficients  are  usually  retained  to  characterize  each 
shape,  the  roots  of  f  (^)  can  be  obtained.  This  will  be  a  very  expensive  com¬ 
putation  however. 

To  reduce  the  expense  in  normalization  several  authors  ([WALL80], 
fKUIIL],  [GRAN72])  have  used  a  "suboptimal"  method  vdiich  normalizes  each 
shape  to  a  standard  orientation  independently.  Because  the  fundamental 


or  c  ]  component  has  been  observed  always  to  be  the  largest  for  a  closed 

curve  traced  CCW,  this  component  is  used  for  scale  normalization.  This  is 

cn 

performed  by  scaling  the  shape  so  that  c',  =  1,  i.e.  c'n  =  p  Most  of  the 

|cil 

variation  between  normalization  methods  comes  in  normalizing  with  respect 
to  rotation  and  starting  point.  The  simplest  method  is  to  rotate  and  shift  so 
that  the  phase  of  the  transformed  coefficients  c  x  and  c  are  zero.  This  is 
equivalent  to  assuming  the  shape  is  an  ellipse  and  rotating  and  shifting  the 
starting  point  so  that  the  major  axis  of  the  ellipse  coincides  with  the  x-axis. 

There  are  two  basic  problems  with  this  method  of  normalization.  For 
many  shapes,  including  all  that  have  N-fold  symmetry  for  N> 2,  we  have  that 
|  c  _|  |  =0,  except  for  representation  errors. 

Secondly,  there  are  two  possible  rotations  of  the  object  that  will  satisfy 
the  above  requirements.  A  180°  rotation  and  T/Z  shift  of  starting  point  also 
places  the  major  axis  of  the  fundamental  ellipse  along  the  x-axis. 

WaUace[WALL80]  addresses  both  of  these  problems.  Wallace’s  normali¬ 
zation  method  is  as  follows: 

1)  Setc0  =  0. 

Z)  Scale  the  object  so  that  |  c ,  j  =  1. 

3)  Find  k  such  that  |  c*  |  is  a  maximum  for  all  k  *0,  1. 

4)  Rotate  the  object  so  that  =  0  and  -  0,  i.e. 

,  _  »[(»-*)* ,  +  (l-n)£*j/(fc-l) 

c  n  —  e 

3)  There  are  m[k]  -1=  |A;  — 1|  —  1,  possible  normalizations  such  that  4)  is 
true.  So,  choose  the  normalization  such  that  £Rejc'nj|Re  [c'n]  I  is  a 

maximum.  This  is  the  normalization  that  places  the  starting  point  and 
the  lobe  of  the  m[k]-fold  fundamental  shape  on  the  x-axis  such  that  the 
curve  is  farthest  from  the  centroid. 


III.  Description  of  System  Implementation 
A.  Library  Generation 

The  first  order  of  business  in  implementing  a  shape  recognition  system 
is  to  generate  the  libraries  for  the  different  shapes.  To  accomplish  this  a 
set  of  parameters  for  an  airplane  are  provided  to  the  program  pLanen2  This 
program  takes  a  parametric  or  "blueprint"  description  of  the  airplane  and 
generates  planar  patches  covering  the  skin  of  the  airplane.  The  parameters 
needed  for  planen2  are  such  things  as:  the  number  of  engines,  where  the 
engines  are  on  the  airplane,  the  length  and  width  of  the  fuselage,  the  posi¬ 
tion  of  the  wings  on  the  fuselage,  the  position  of  the  wing  tip  edges,  the 
width  of  the  wings  at  the  fuselage  and  at  the  wings  tips,  the  same  data  for 
the  tail,  wings,  and  tail  is  needed.  The  fuselage  and  engine  bodies  arc 
approximated  by  a  cylindrical  prism.  The  nose  cone  and  tail  cone  on  the 
fuselage  are  approximated  by  a  prismatic  cone.  Coordinates  for  the  ver¬ 
tices  of  the  planar  patches  are  output  fro'm  the  program. 

For  each  aspect  desired  in  the  library  for  each  airplane,  the  program 
mdgen  generates  a  chain-code.  The  chain-code  for  the  desired  aspect  angle 
is  obtain  as  follows: 

a)  Rotate  each  vertex  of  each  planar  patch  by  the  aspect  angle  ,  <pyg) 

b)  Project  the  rotated  patch  onto  the  image  plane  (set  z  coordinate  to 
zero). 

c)  The  object  is  traced  and  the  chain-code  is  generated  by  detecting  if  a 
point  in  a  NxN  grid  in  the  image  plane  lies  within  at  least  one  of  the  pro¬ 
jected  planar  patches.  This  part  of  the  algorithm  is  sped  up  by  sorting 
the  planar  patches  by  the  upper  left  hand  corner  coordinate  of  the 
patch  and  its  lower  right  hand  corner  coordinate.  An  initial  edge  point 
on  the  contour  is  found  by  beginning  at  the  center  of  the  grid  and  scan¬ 
ning  to  find  an  object  edge  transition.  Once  a  starting-point  on  the 
edge  is  found,  then  the  nearest  patches  in  the  image  plane  are 
searched  first  to  indicate  if  that  grid  point  lies  within  one  of  the 
patches.  From  this  edge  following,  the  object  contour  is  traced  and  its 
corresponding  chain-code  is  produced. 


Once  chain-code  descriptions  have  been  determined  for  all  the  desired 
aspects  for  each  airplane,  the  normalized  Fourier  descriptor  is  calculated. 


B.  Chain-Code  Format 

The  chain-code  output  from  mdgan  has  with  it  attributes  that  mi#it  be 
used  later  in  processing  the  chain-code.  Because  of  this  a  standard  chain- 
code  format  was  adopted  having  the  following  form: 


short  int 

id 

unique  chain-code  id  number 

short  int 

hags 

hags 

short  int 

code 

image  code 

short  int 

xstart 

x  starting  coordinate 

short  int 

ystart 

y  starting  coordinate 

short  int 

zstart 

z  starting  coordinate 

short  int 

scale 

object  scale 

short int 

rotang 

angle  of  rotation 

short  int 

length 

length  of  chain-code 

char 

chain[length] 

chain-code 

The  chain-code  itself  is  an  eight-neighbor  chain-code  which  indicates  the 
direction  to  proceed  to  the  next  point  along  the  contour.  The  chain-code 
directions  are  encoded  as  follows: 


3  2  1 

4  +  0 

5  6  7 


C.  Normalised  Fourier  Descriptor  Calculation 

The  normalized  Fourier  descriptor  (NFD)  of  each  chain-code  is  calcu¬ 
lated  (is  calculated  using  the  program  mrbLd)  .The  normalized  Fourier 
descriptor  (NFD)  of  each  chain-code  is  calculated  by  carrying  out  the  follow¬ 
ing  operations: 

a)  The  chain-code  is  converted  to  a  corresponding  complex  vector 

f  =  (*y  +  t  yj). 

b)  The  complex  vector  is  smoothed  by  using  a  convolutional  window  as  a 
function  of  arc-length.  Possible  windows  are  rectangular  ,  triangular 
and  Gaussian. 

c)  Next,  if  necessary,  the  contour  complex  vector  is  resampled  to  a  power 
of  two  using  linear  interpolation  (with  respect  to  arc-length). 

d)  Once  we  have  a  complex  vector,that  is  a  power  of  two,  we  can  perform  a 
radix  two  FFT. 

e)  The  result  of  the  FFT  is  then  normalized  using  the  sub-optimal  method 
of  Wallace  described  above.  Both  the  NFD  ,  the  starting-point  rotation 
angle,  and  a  normalization  signature  are  output  by  this  program.  The 
normalization  signature  provides  the  information  concerning  the 
coefficient  used  in  the  normalization  and  which  lobe  was  rotated  to  the 
x-axis. 


D.  Feature  Vector  Format 

Later,  the  output  of  the  NFD  generation  program  will  be  in  a  standard 
"feature  vector  "  format.  A  standard  feature  format  is  desired  so  that  com¬ 
peting  shape  methods  can  be  analyzed  using  the  same  feature  reduction 
and  classification  programs.  This  format  is  essentially  the  same  as  the 
chain-code  format  except  in  the  PDS  header  information  is  provided  to 
designate  what  method  was  used  to  create  this  feature  vector  and  ■what  is 
the  arithmetic  type  of  the  data  vector  (i.e.  character,  integer,  floating,  or 
complex  data).  At  the  present,  the  NFDs  is  written  in  Fortran  formatted 
ASCII  I/O. 


E.  Feature  Reduction 

The  features  provided  by  the  different  shape  methods  is  assumed  to 
contain  some  redundant  information,  so  some  kind  of  feature  vector 
transformation  is  desired  to  reduce  the  dimensionality  of  the  feature  vec¬ 
tor.  This  is  necessary  to  reduce  the  storage  requirements  for  the  libraries . 
This  reduction  is  accomplished  by  performing  an  eigen-analysis  (using  the 
program  bldeig )  on  all  the  libraries.  The  feature  vectors  are  eigen- analyzed 
by  forming  the  matrix  A(i,j)  - 1  tT  and  averaging  over  all  features  in  all 
the  libraries.  The  resulting  eigenvectors  are  sorted  by  there  eigenvalues. 
The  eigenvectors  corresponding  to  the  K  largest  eigenvalues  are  used  to 
transform  feature  vectors. 

The  feature  vectors  in  each  library  are  transformed  by  the  eigenvectors 
having  the  largest  corresponding  eigenvalues  as  described  above  using  the 
program  bldint.  Each  component  is  multiplied  by  10,000  and  the  integer 
part  retained  in  order  to  reduce  the  storage  requirement  for  the  libraries . 
The  original  feature  vectors  were  Nx 4  bytes  and  the  transformed  feature 
vectors  are  Kx2  bytes  (K<N). 

F.  Unknown  Feature  Vector  Extraction 

To  obtain  the  unknown  feature  vector,  first  the  object  must  be 
extracted  or  segmented  from  the  background.  This  is  a  very  image  class 
dependent  problem  and  many  methods  exist  to  accomplish  this.  For  the 
pan-chromatic  video  data  provided,  in  order  to  segment  the  object  from  the 
background,  a  rather  sophisticated  image  segmentation  algorithm  using 
both  the  grey-levels  and  edge  infonr  tion  was  required.  These  operations 
are  performed  by  the  programs  to  kef,  pottarg2,  mean,  mountain,  trimp- 
targ2,  and  radts*g2. 

Once,  the  object  is  segmented  from  the  background,  the  object  boun¬ 
dary  is  traced  and  its  corresponding  chain-code  output.  The  tracing  is  per¬ 
formed  by  the  program  ctroce.  The  chain-code  is  input  to  the  same  pro¬ 
grams  that  generated  the  feature  vectors  for  the  library,  except  in  this  case 
it  is  often  necessary  to  filter  the  contour  more  for  noise  removal. 
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G.  Feature  Vector  Classification. 

Once  the  feature  vector  (i.e.  the  transformed  NFD)  has  been  computed, 
this  unknown  feature  vector  is  compared  to  all  the  entries  in  each  library. 
The  unknown  feature  vector  is  then  assigned  to  the  airplane  and  to  that 
aspect  for  which  the  distance  between  the  unknown  feature  vector  and  the 
library  feature  vector  is  the  smallest.  In  our  classification  program,  the 
sum  of  the  squared  differences  between  the  components  of  the  feature  vec¬ 
tors  is  used  as  the  distance  measure.  This  distance  measure  is  used  since  it 
corresponds  to  the  sum-of-squares  distance  measure  in  the  space  domain 
as  well.  The  classification  is  performed  by  the  program  clasfy. 


IV.  Results  of  Fourier  Descriptors  Applied  to  Target  Recognition 

The  algorithms  described  above  were  used  to  perform  recognition 
experiments  on  the  data  supplied  by  RADC.  (This  data  is  described  in  the 
final  report  for  SCEE  Contract  PDP81-31.)  The  results  indicate  that  when  the 
target  is  sufficiently  different  in  gray  level  from  the  background,  the  shape 
can  be  accurately  extracted  and  very  high  (>97%)  classification  accuracies 
can  be  obtained.  However  under  more  commmon  conditions,  the  lighting, 
shadows,  and  background  variations  cause  difficulty  with  object  segmenta¬ 
tion,  and  the  classifier  must  be  adjusted  to  allow  more  variation  in  the 
expected  contour  (we  allow  12%  mean-square-error).  In  this  case, 
classification  accuracy  drops  to  90%  for  a  four  class  problem  and  false 
alarms  occur  approximately  once  per  512x512  pixel  image. 
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